***load data

use "${path}\data\naep\naep_ready.dta", clear 

***Main Tables

cd "${path}\output\tables\"

***Table I
*Panel A (Run GSS and ACS Analysis afterwards for Panels B and C)

eststo: reg evo_correct evo_score if schtype==1 & sample_naep==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls  if schtype==1  & sample_naep==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score i.state i.byear if schtype==1  & sample_naep==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1  & sample_naep==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
esttab using main_table_all.tex, replace keep(evo_score) b(3) se(3) fragment  nomtitles nonotes  label ///
stats(ymean ysd  r2_a NObs, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  /// 
star(* 0.1 ** 0.05 *** 0.01)  ///
refcat(evo_score "\specialcellleft{\textbf{\emph{Panel A: Outcome: Evolution Knowledge in School}}}", nolabel) ///
mlabels("\specialcell{Controls : NO \\ State FE: NO \\ Cohort FE: NO}" "\specialcell{Controls : YES \\ State FE: NO \\ Cohort FE: NO}" "\specialcell{Controls : NO \\ State FE: YES \\ Cohort FE: YES}" "\specialcell{Controls : YES \\ State FE: YES \\ Cohort FE: YES}") 
eststo clear

***Table II
*Panel A (Run GSS and ACS Analysis afterwards for Panels B and C)

eststo: reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1 & closestate10==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls democrat i.state i.year i.byear if schtype==1 , vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls_trend $statelineartrend $sq_statelineartrend i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
rename evo_score evo_score_temp
rename evo_score_dosage evo_score
eststo: reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
rename evo_score evo_score_dosage
rename evo_score_temp evo_score
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & onereform==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & largestate==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & textstate==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & highschoolentryyear>1994, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & highschoolentryyear>1999, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
probit evo_correct evo_score $controls i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
estadd scalar r2_a =round(e(r2_p), .001)
estpost margins, dydx(evo_score)
eststo ames2
eststo: reg v2_evo_correct evo_score $controls i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
esttab using further_robustness_all.tex, replace keep(evo_score) b(3) se(3) fragment noobs nomtitles nonotes  label ///
mlabels("\specialcell{Close \\ Elections}" "\specialcell{Control: \\ Governor's \\ Party}" "\specialcell{State Specific \\ Time Trends}" "\specialcell{Dosage \\ Treatment}" "\specialcell{Only One \\ Reform Event}" "\specialcell{Only Large \\ States}" "\specialcell{Only States \\ With Std. Text}""\specialcell{Sample Start: \\ 1995}" "\specialcell{Sample Start: \\ 2000}" "Probit" "\specialcell{Outcome Coding: \\ Variation 1}" "\specialcell{Outcome Coding \\ Variation 2}") ///
refcat(evo_score "\specialcellleft{\textbf{\emph{Panel A: Outcome: Evolution Knowledge in School}}}", nolabel) ///
star(* 0.1 ** 0.05 *** 0.01)  
eststo clear

***Online Appendix Tables

cd "${path}\output\tables\"

***Table A.III

eststo: reg private evo_score $controls_pub_priv i.state i.year i.byear if sample_naep_priv_pub==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg private evo_score $controls_pub_priv i.state i.year i.byear if female==1 & sample_naep_priv_pub==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg private evo_score $controls_pub_priv i.state i.year i.byear if pc_athome==0 & sample_naep_priv_pub==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls_pub_priv i.state i.year i.byear if sample_naep_priv_pub==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls_pub_priv  i.state i.year i.byear if schtype==1 & sample_naep_priv_pub==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls_pub_priv i.state i.year i.byear if schtype>1 & sample_naep_priv_pub==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
esttab using naep_evolution_by_schooltype.tex, replace keep(evo_score) b(3) se(3)  fragment  nonotes  label ///
stats(ymean ysd  r2_a NObs, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  /// 
mlabels("\specialcell{All \\ Students}" "\specialcell{Only \\ Female}" "\specialcell{Only Without \\ PC At Home}" "\specialcell{All \\ Students}" "\specialcell{Only Public \\ School Students}" "\specialcell{Only Private \\ School Students}") ///
mgroups("Attending Private School" "Evolution Knowledge", pattern(1 0 0 1 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) /// 
star(* 0.1 ** 0.05 *** 0.01)  ///
indicate("State FEs = *state" "Cohort FEs = *byear" "Controls = female", labels("YES" "NO"))
eststo clear

***Table A.IV
*Panel A

centile (evo_score) if sample_less==1, centile (20 80)

csdid evo_correct_conv if sample_less==1, time(highschoolentryyearcs) gvar(first_treat) notyet agg(event) method(dripw) rseed(2) cluster(state)
estat simple

csdid evo_correct_conv if sample_less==1 & evo_score<0.91 & evo_score>0.64, time(highschoolentryyearcs) notyet gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(state)
estat simple

centile (evo_score) if sample_more==1, centile (20 80)

csdid evo_correct_conv if sample_more==1, time(highschoolentryyearcs) gvar(first_treat) notyet agg(event) method(dripw) rseed(2) cluster(state)
estat simple

csdid evo_correct_conv if sample_more==1 & evo_score<0.82 & evo_score>0.07, time(highschoolentryyearcs) notyet gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(state)
estat simple


***Table A.V

eststo: reg biology1_ever evo_score $controls  i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg biology2_ever evo_score $controls  i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg biologyap_ever evo_score $controls  i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
esttab using naep_mechanism_course_selection.tex, replace keep(evo_score) b(3) se(3) fragment  nonotes  label /// 
stats(ymean ysd  r2_a NObs, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  /// 
mlabels("\specialcell{First Year \\ Biology}" "\specialcell{Second Year \\ Biology}" "\specialcell{Advanced \\ Placement \\ Biology}")  ///
star(* 0.1 ** 0.05 *** 0.01)  ///
indicate("State FEs = *state" "Cohort FEs = *byear" "Controls = female", labels("YES" "NO"))
eststo clear

***Table A.VI

foreach v of varlist  evo_correct placebo_share_correct_all motion_correct matterandmass_correct energy_correct reproduction_correct climate_correct pollution_correct earth_correct tectonics_correct universe_correct  {
eststo: reg `v' evo_score $controls  i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
}
esttab using naep_placebo_knowledge.tex, replace keep(evo_score) b(3) se(3) fragment  nonotes  label /// 
stats(ymean ysd  r2_a NObs, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  /// 
mgroups("Main Outcome:" "Non-Evolution Scientific Topics:", pattern(1 1 0 0 0 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) /// 
mlabels("Evolution" "\specialcell{Average}" "Motion" "\specialcell{Matter \\ and Mass}" "Energy" "Reproduction" "\specialcell{Climate \\ Change}" "Pollution" "Earth" "Tectonics" "Universe") /// 
star(* 0.1 ** 0.05 *** 0.01)  
eststo clear

***Table A.XII 

eststo: reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1 & closestate10==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls democrat i.state i.year i.byear if schtype==1 , vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls_trend $statelineartrend $sq_statelineartrend i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
rename evo_score evo_score_temp
rename evo_score_dosage evo_score
eststo: reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
rename evo_score evo_score_dosage
rename evo_score_temp evo_score
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & onereform==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & largestate==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & textstate==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & highschoolentryyear>1994, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
eststo: reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & highschoolentryyear>1999, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
probit evo_correct evo_score $controls i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
estadd scalar r2_a =round(e(r2_p), .001)
estpost margins, dydx(evo_score)
eststo ames2
eststo: reg v2_evo_correct evo_score $controls i.state i.year i.byear if schtype==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
esttab using naep_further_robustness.tex, replace keep(evo_score) b(3) se(3) fragment nomtitles nonotes  label ///
stats(ymean ysd r2_a NObs, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  ///
mlabels("\specialcell{Close \\ Elections}" "\specialcell{Control: \\ Governor's \\ Party}" "\specialcell{State Specific \\ Time Trends}" "\specialcell{Dosage \\ Treatment}" "\specialcell{Only One \\ Reform Event}" "\specialcell{Only Large \\ States}" "\specialcell{Only States \\ With Std. Text}""\specialcell{Sample Start: \\ 1995}" "\specialcell{Sample Start: \\ 2000}" "Probit" "\specialcell{Outcome Coding: \\ Indicator Variation}") ///
star(* 0.1 ** 0.05 *** 0.01)  ///
indicate("State FEs = *state" "Cohort FEs = *byear" "Controls = female" , labels("YES" "YES"))
eststo clear

***test for statistical significance

fvset base 5 state
fvset base 1990 byear
fvset base 2009 year

reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1
estimates store base

reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1 & closestate10==1
estimates store close

reg evo_correct evo_score $controls democrat i.state i.year i.byear if schtype==1
estimates store demo

reg evo_correct evo_score $controls $statelineartrend $sq_statelineartrend i.state i.year i.byear if schtype==1
estimates store trends

rename evo_score evo_score_temp
rename evo_score_dosage evo_score
reg evo_correct evo_score $controls  i.state i.year i.byear if schtype==1
estimates store dosage
rename evo_score evo_score_dosage
rename evo_score_temp evo_score

reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & onereform==1
estimates store oneref

reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & largestate==1
estimates store largestate

reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & textstate==1
estimates store textstate

reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & highschoolentryyear>1994
estimates store samp1

reg evo_correct evo_score $controls i.state i.year i.byear if schtype==1 & highschoolentryyear>1999
estimates store samp2

reg v2_evo_correct evo_score $controls i.state i.year i.byear if schtype==1
estimates store outc1

suest base close, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[close_mean:evo_score]

suest base demo, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[demo_mean:evo_score]

suest base trends, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[trends_mean:evo_score]

suest base dosage, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[dosage_mean:evo_score]

suest base oneref, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[oneref_mean:evo_score]

suest base largestate, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[largestate_mean:evo_score]

suest base textstate, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[textstate_mean:evo_score]

suest base samp1, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[samp1_mean:evo_score]

suest base samp2, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[samp2_mean:evo_score]

suest base outc1, vce(cluster state16) coeflegend
test _b[base_mean:evo_score] = _b[outc1_mean:evo_score]

quietly generate yreg = evo_correct
quietly generate yprobit = evo_correct
gsem (yprobit  <-  evo_score $controls  i.state i.year i.byear, probit  vce(cluster state16)) (yreg <-  evo_score $controls  i.state i.year i.byear if schtype==1)
margins, dydx(evo_score) post
margins, coeflegend
test  _b[evo_score:1bn._predict]= _b[evo_score:2._predict]

***Table A.XIX

preserve
	label variable evo_correct "Evolution Knowledge"
eststo: estpost corr evo_correct placebo_share_correct_all motion_correct matterandmass_correct energy_correct reproduction_correct climate_correct pollution_correct earth_correct tectonics_correct universe_correct  
esttab using naep_scientific_areas_correlation_short.tex, nonumber not unstack compress noobs replace fragment nonotes label star(* 0.1 ** 0.05 *** 0.01)  
eststo clear
restore

***Table A.XX

preserve
foreach v of varlist * {
	label variable `v' `"\hspace{0.1cm} `: variable label `v''"'
	}
eststo: estpost sum evo_score evo_correct placebo_share_correct_all motion_correct matterandmass_correct energy_correct reproduction_correct climate_correct pollution_correct earth_correct tectonics_correct universe_correct female race_white race_black race_hispanic race_asian race_other sublunch books_athome_few books_athome_shelf books_athome_bookcase books_athome_many pc_athome sessno bmonth if schtype==1 
esttab using naep_long_descriptive_statistics.tex, replace wide cells("mean(fmt(%9.2fc) label(Mean)) sd(fmt(%9.2fc) label(\specialcell{Std. \\ Dev.})) min(fmt(%9.2fc) label(Min.)) max(fmt(%9.2fc) label(Max.))")  label booktabs nonum fragment noobs  ///
refcat(evo_score "\emph{Treatment Variable:}" evo_correct "\emph{Main Outcome:}" placebo_share_correct_all "\emph{Non-Evolution Scientific Outcomes:}" female "\emph{Controls:}", nolabel)
eststo clear
restore

***Table A.XXXII

foreach i of num 9/1 { 
eststo: reg evo_correct evo_score_`i' $controls  i.state i.year i.byear if schtype==1 & sample_naep==1, vce(cluster state) baselevels
estadd ysumm
estadd scalar NObs =round(e(N),10)
}
esttab using naep_by_level.tex, replace keep(evo_score_9 evo_score_8 evo_score_7 evo_score_6 evo_score_5 evo_score_4 evo_score_3 evo_score_2 evo_score_1) b(3) se(3) fragment nomtitles nonotes  label ///
stats(ymean ysd r2_a NObs, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations")  fmt(2 2 3 %9.0fc))  ///
mgroups("Evolution Knowledge", pattern(1 0 0 0 0 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) /// 
indicate("State FEs = *state" "Cohort FEs = *byear" "Controls = female" , labels("YES" "NO")) ///
star(* 0.1 ** 0.05 *** 0.01)  
eststo clear   

***Main Figures
cd "${path}\output\figures"

***Figure II
*Panel A

reg evo_correct_conv pre78 pre56 pre34 zero post01 post23 post45 $controls i.year i.state i.byear if sample_less==1, vce(cluster state)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin)) alpha(0.05)

graph save "event_study_TWFE.pdf", replace

csdid evo_correct_conv if sample_less==1, time(highschoolentryyearcs) gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(state)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))  alpha(0.05)

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together  noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Evolution Knowledge") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\naep-less-twfecs.pdf", replace

***Figure II
*Panel B

reg evo_correct_conv pre78 pre56 pre34 zero post01 post23 post45 $controls i.year i.state i.byear if sample_more==1, vce(cluster state)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin)) alpha(0.05)

graph save "event_study_TWFE.pdf", replace

csdid evo_correct_conv if sample_more==1, time(highschoolentryyearcs) gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(state)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))  alpha(0.05)

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together  noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Evolution Knowledge") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\naep-more-twfecs.pdf", replace


***Online Appendix Figures

cd "${path}\output\figures"

***Figure A.III
*Panel B

foreach v of varlist sublunch books_athome_few books_athome_shelf bookcase books_athome_many pc_athome sessno bmonth {
eststo: reg evo_score `v' i.year if sample_naep==1, cluster(state16) 
estimates store dose_`v'
}

coefplot (dose_sublunch, keep(sublunch)) (dose_books_athome_few, keep(books_athome_few)) (dose_books_athome_shelf, keep(books_athome_shelf)) (dose_bookcase, keep(bookcase)) (dose_books_athome_many, keep(books_athome_many)) (dose_pc_athome, keep(pc_athome)) (dose_sessno, keep(sessno)) (dose_bmonth, keep(bmonth)), nokey graphregion(color(white)) bgcolor(white) xline(0) xlabel(-0.3(0.1)0.3) ciopts(color(black)) mcolor(black) msymbol(o) msize(small)
graph export "${path}\output\figures\dose_naep.pdf", replace
eststo clear

***Figure A.IV

preserve
cd "${path}\data\naep\"
clear
set obs 38
set seed 2
 
 foreach num of numlist 1/1000 {
  gen reformyear`num'=.
  replace reformyear`num' = 2005 in 1
  replace reformyear`num' = 2006 in 2
  replace reformyear`num' = 2005 in 3
  replace reformyear`num' = 2009 in 4
  replace reformyear`num' = 2004 in 5
  replace reformyear`num' = 2006 in 6
  replace reformyear`num' = 2006 in 7
  replace reformyear`num' = 2008 in 8
  replace reformyear`num' = 2004 in 9
  replace reformyear`num' = 2005 in 10
  replace reformyear`num' = 2004 in 11
  replace reformyear`num' = 2006 in 12
  replace reformyear`num' = 2007 in 13
  replace reformyear`num' = 2007 in 14
  replace reformyear`num' = 2005 in 15
  replace reformyear`num' = 2007 in 16
  replace reformyear`num' = 2002 in 17
  replace reformyear`num' = 2006 in 18
  replace reformyear`num' = 2000 in 19
  replace reformyear`num' = 2009 in 20
  replace reformyear`num' = 2008 in 21
  replace reformyear`num' = 2008 in 22
  replace reformyear`num' = 2006 in 23
  replace reformyear`num' = 2004 in 24
  replace reformyear`num' = 2006 in 25
  replace reformyear`num' = 2003 in 26
  replace reformyear`num' = 2004 in 27
  replace reformyear`num' = 2006 in 28
  replace reformyear`num' = 2006 in 29  
  replace reformyear`num' = 2002 in 30
  replace reformyear`num' = 2006 in 31
  replace reformyear`num' = 2005 in 32
  replace reformyear`num' = 2005 in 33
  replace reformyear`num' = 2007 in 34
  replace reformyear`num' = 2009 in 35
  replace reformyear`num' = 2003 in 36
  replace reformyear`num' = 2008 in 37
  replace reformyear`num' = 2003 in 38
 
  
   gen rand`num' = runiform()
  sort rand`num'
  drop rand`num'
 
 }

 gen state=.
 replace state= 1 in 1
 replace state= 2 in 2
 replace state= 5 in 3
 replace state= 8 in 4
 replace state= 9 in 5
 replace state= 10 in 6
 replace state= 11 in 7
 replace state= 12 in 8
 replace state= 13 in 9
 replace state= 15 in 10 
 replace state= 17 in 11
 replace state= 18 in 12
 replace state= 19 in 13
 replace state= 20 in 14
 replace state= 22 in 15
 replace state= 23 in 16
 replace state= 24 in 17
 replace state= 25 in 18
 replace state= 26 in 19
 replace state= 27 in 20
 replace state= 28 in 21
 replace state= 29 in 22
 replace state= 30 in 23
 replace state= 32 in 24
 replace state= 33 in 25
 replace state= 35 in 26
 replace state= 37 in 27
 replace state= 38 in 28
 replace state= 39 in 29
 replace state= 42 in 30
 replace state= 44 in 31 
 replace state= 45 in 32
 replace state= 46 in 33
 replace state= 47 in 34
 replace state= 48 in 35
 replace state= 51 in 36
 replace state= 54 in 37
 replace state= 56 in 38
 
order state
save "permutation.dta", replace 
restore

merge m:1 state using "permutation.dta" 

tempname myresults 
postfile `myresults' threshold intercept gradient se using myresults.dta, replace 
 foreach i of num 1/1000 {
 
gen evolution_score_perm_`i'=.
replace evolution_score_perm_`i'=9 if state==1 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=21 if state==1 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=48 if state==2 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=59 if state==2 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==4 
replace evolution_score_perm_`i'=55 if state==5 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=66 if state==5 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==6
replace evolution_score_perm_`i'=86 if state==8 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==8 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==9 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=59 if state==9 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=91 if state==10 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=80 if state==10 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=80 if state==11 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==11 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=16 if state==12 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==12 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=7 if state==13 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=66 if state==13 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=91 if state==15 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=75 if state==15 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==16 
replace evolution_score_perm_`i'=45 if state==17 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==17 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==18 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==18 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=77 if state==19 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=0 if state==20 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==20 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=55 if state==21 
replace evolution_score_perm_`i'=64 if state==22 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=27 if state==22 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=30 if state==23 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=68 if state==23 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=77 if state==24 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=73 if state==24 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==25 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=84 if state==25 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=84 if state==26 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=80 if state==26 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=86 if state==27 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=89 if state==27 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=5 if state==28 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=87 if state==28 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==29 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=78 if state==29 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==30 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=75 if state==30 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=66 if state==31
replace evolution_score_perm_`i'=70 if state==32 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=77 if state==32 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=23 if state==33 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==33 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==34 
replace evolution_score_perm_`i'=73 if state==35 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==35 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=68 if state==36
replace evolution_score_perm_`i'=100 if state==37 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==37 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=9 if state==38 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=64 if state==38 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=30 if state==39 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=86 if state==39 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=25 if state==40 
replace evolution_score_perm_`i'=82 if state==41 
replace evolution_score_perm_`i'=91 if state==42 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==42 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==44 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==44 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=95 if state==45 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==45 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==46 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=77 if state==46 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=2 if state==47 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=55 if state==47 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=64 if state==48 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=46 if state==48 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==49 
replace evolution_score_perm_`i'=86 if state==50 
replace evolution_score_perm_`i'=50 if state==51 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=68 if state==51 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=86 if state==53 
replace evolution_score_perm_`i'=2 if state==54 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=46 if state==54 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=55 if state==55
replace evolution_score_perm_`i'=36 if state==56 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=61 if state==56 & highschoolentryyear>=reformyear`i'
gen evo_score_perm_`i' = evolution_score_perm_`i'/100
label variable evo_score_perm_`i' "Evolution Score"

reg evo_correct evo_score_perm_`i' $controls  i.state i.year i.byear if sample_naep==1 & schtype==1, vce(cluster state) baselevels
post `myresults' (`i') (`=_b[_cons]') (`=_b[evo_score_perm_`i']') (`=_se[evo_score_perm_`i']') 
}
postclose `myresults' 

preserve
use "${path}\data\naep\myresults.dta", clear
twoway kdensity gradient, lcolor(black) xline(0.0648, lcolor(red)) xline(0.0606, lcolor(black) lpattern(dash)) graphregion(fcolor(white)) bgcolor(white) xtitle("Coefficients from placebo regressions") ytitle("Density")
graph export "${path}\output\figures\permutations_naep.pdf", replace
sum gradient, d
restore

***Figure A.VII
*Panel A

reghdfe evo_correct_conv pre78 pre56 pre34 zero post01 post23 post45 $controls if sample_less==1, absorb(year state byear) cluster(state) noconstant
matrix l_vec = 1 \ 0 \ 0
local plotopts xtitle(Mbar) ytitle(95% Robust CI)
honestdid, l_vec(l_vec) pre(1/3) post(4/6) mvec(0(0.0005)0.004) omit alpha(0.05) delta(sd) coefplot `plotopts' yscale(r(0 0.1)) ytick(0(0.02)0.1) ylabel(0(0.02)0.1)
graph export "${path}\output\figures\naep-less-rr.pdf", replace

***Figure A.VII
*Panel B

reghdfe evo_correct_conv pre78 pre56 pre34 zero post01 post23 post45 $controls if sample_more==1, absorb(year state byear) cluster(state) noconstant
matrix l_vec = 1 \ 0 \ 0
local plotopts xtitle(Mbar) ytitle(95% Robust CI)
honestdid, l_vec(l_vec) pre(1/3) post(4/6) mvec(0(0.0005)0.004) omit alpha(0.05) delta(sd) coefplot `plotopts' yscale(r(-0.02 0.08)) ytick(-0.02(0.02)0.08) ylabel(-0.02(0.02)0.08)
graph export "${path}\output\figures\naep-more-rr.pdf", replace

***Figure A.X
*Panel A

reg evo_correct_conv pre78 pre56 pre34 zero post01 post23 post45 $controls i.year i.state i.byear if sample_less==1 | state==31 | state==41 | state==49 | state==55, vce(cluster state)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin)) alpha(0.1)

graph save "event_study_TWFE.pdf", replace

csdid evo_correct_conv if sample_less==1  | state==31 | state==41 | state==49 | state==55, time(highschoolentryyearcs) gvar(first_treat) notyet agg(event) method(dripw) rseed(2) cluster(state)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))  alpha(0.1)

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together  noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Evolution Knowledge") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\naep-less-twfecs_nevert.pdf", replace

***Figure A.X
*Panel B

reg evo_correct_conv pre78 pre56 pre34 zero post01 post23 post45 $controls i.year i.state i.byear if sample_more==1 | state==31 | state==41 | state==49 | state==55, vce(cluster state)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin)) alpha(0.1)

graph save "event_study_TWFE.pdf", replace

csdid evo_correct_conv if sample_more==1  | state==31 | state==41 | state==49 | state==55, time(highschoolentryyearcs) gvar(first_treat) notyet agg(event) method(dripw) rseed(2) cluster(state)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(-0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.2 0.2)) ytick(-0.2(0.1)0.2) ///
		ylabel(-0.2(0.1)0.2, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Evolution Knowledge", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))  alpha(0.1)

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together  noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Evolution Knowledge") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\naep-more-twfecs_nevert.pdf", replace

***Figure A.XIV

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1, vce(cluster state) baselevels
estimates store Overall

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & female==1, vce(cluster state) baselevels
estimates store Female

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & female==0, vce(cluster state) baselevels
estimates store Male

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & race_white==1, vce(cluster state) baselevels
estimates store White

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & race_black==1, vce(cluster state) baselevels
estimates store Black

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & race_hispanic==1, vce(cluster state) baselevels
estimates store Hispanic

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & sublunch==1, vce(cluster state) baselevels
estimates store Subsidized

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & sublunch==0, vce(cluster state) baselevels
estimates store Non_Subsidized

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & pc_athome==1, vce(cluster state) baselevels
estimates store PC

reg evo_correct evo_score $controls i.state i.year i.byear if sample_naep==1 & pc_athome==0, vce(cluster state) baselevels
estimates store No_PC

coefplot Overall || Female || Male || White || Black || Hispanic || Subsidized || Non_Subsidized || PC || No_PC, keep(evo_score)  xline(0) bycoefs byopts(xrescale) graphregion(color(white)) bgcolor(white) levels(95) bylabels("Overall (n=14,080)" "Female (n=7,280)" "Male (n=6,800)" "White (n=7,990)" "Black (n=2,650)" "Hispanic (n=2,390)" "Subsidized (n=4,270)" "Not Subsidized (n=9,810)" "PC (n=12,480)" "No PC (n=1,600)") headings(2 = "{bf:By Gender:}" 4= "{bf:By Race/Ethnicity:}"  7= "{bf:By Subsidized Lunch Status:}" 9= "{bf:By PC at Home:}") 
    graph export "${path}\output\figures\naep_coefplot_heterogeneities.pdf", replace